//To compile in PuTTy, have this file, and VerificationSim.h in the same folder, and
//type: "g++ VerificationSim.cpp".

#include "VerificationSim.h"

//#include <iostream>

//Verification Sim
using namespace std;

Matrix * A;
Matrix * B;
Matrix * C;
//extern int n;

bool trace;

int main(){
	time_t seconds;

	seconds = time (NULL);
//Variables
   double lambda = 15;
   lambda = lambda /14;
   Qbit z;
   bool equal = true;
   srand(time(0)); // Seed the random number generator
   Matrix pauliZ = Matrix::SetPauliZ();
   int n;
   

//Introduction, and Inputs   
   cout << "Algorithm Simulator for the Verification of Matrix products" << endl;
   cout << "Please input n, the size of the matrices: ";
   cin >> n;
   double start = 2;
	start = start/3;
   double looplimit = (log(pow((double)n,start))/log(lambda)) +9; //using change-of-base formula
   //cout << "LOOPLIMIT" << looplimit;
   cout << "Please input the values of the matrix A (L to R by row)" << endl;
   //initialize A
   A = new Matrix(n, true);
   
   cout << endl << "Please now input the values of the matrix B";
   //initialize B
   B = new Matrix(n, true);
   
   cout << endl << "Please now input the values of the matrix C";
   //initialize C
   C = new Matrix(n, true);
   
   
   //Determine whether the user wants to track the algorithm or not.
   while(true){
      cout << "Would you like to trace the algorithm, and see it step-by-step  <y/n>";
      string a;
      cin >> a;
      if(strcmp(a.c_str(),"y") == 0){
         trace = true;
         break;
      }else
      if(strcmp(a.c_str(),"n")==0){
         trace = false;
         break;
      }
      else
      cout << "--Invalid input--" << endl;
   }


//Optional reprinting of matrices
if(trace){
   A->Print(); cout << endl;
   B->Print(); cout << endl;
   C->Print(); cout << endl;
  }
  
//Now, we run the main loop:
if(trace){ cout << "We now start the algorithm loop running it " << looplimit << " times" << endl;}

//start timer
 clock_t begin=clock();
 
 for(int i = 0; i < looplimit; i++){
      for(int j = 0; j < 16; j++){
         int k = (int)(2*pow(lambda, i));
         int l = rand() % (k+1);
         
		 z = initialization(A, B, C, n, l, k);
		 
		 z.ApplyHgate();
         if(trace){ cout << "We apply a hadamard gate to the result of VerifyOnce()" << endl;}
         z.Normalize();
         if(trace){ cout << "We normalize the result" << endl;}
         z.Measure();
         if(trace){ cout << "And then we measure. " << endl << " If the result is 1, the matrices are not equal and the loop stops" << endl << "if the result is 0, the algorithm loop continues";}

         if(z.re1 == 1){
            equal = false;
           if(trace){ cout << "the result is 1, so we break the loop" << endl;
		   }
		   break;
         }
         if(trace){ cout << "the result is 0, so we continue the loop" << endl;}
      }
	  if(z.re1 == 1){
            equal = false;
		   break;
         }
   }
   
   //stop timer
   clock_t end=clock();
   cout << endl << "Time elapsed: " << (end-start)/double(CLOCKS_PER_SEC)*1000 << " ms"<< endl;
//Output result
   if(equal){
      cout << "The matrices are equal";
   }
   else{
      cout << "The matrices are NOT equal";
   }
 
   return 0;
}


Qbit initialization(Matrix * A, Matrix * B, Matrix * C, int n, int l, int k){
   
   vector<double>p(n);
   vector<double>q(n);
   // we manually set the size of subset to be 3
   int p1;
   int p2;
   int p3; 

   //choose the three vectors in the subset of p
   do{
      p1 = rand() % n;
      p2 = rand() % n;
      p3 = rand() % n;
   }while(! ((p1 != p2) && (p1 != p3) && (p2 != p3))); 
   
   int q1;
   int q2;
   int q3;
   //choose the three vectors in the subset of p
   do{
      q1 = rand() % n;
      q2 = rand() % n;
      q3 = rand() % n;
   }while (!((q1 != q2) && (q1 != q3) && (q2 != q3)));

   //Initialize p and q
   for (int i = 0; i < n; i++){
      p[i] = (double)rand()/RAND_MAX;
      q[i] = (double)rand()/RAND_MAX;
   }
   
   //calculate vector aR, bS,cRS
   vector<double> aR(n);
   vector<double> bS(n);
   vector<double> cr(n);
   vector<double> cs(n);
   double cRS = 0;
   
   for (int j = 0; j < n; j++){		//Segfault in here
       aR[j] = A->Get(p1,j) * p[p1] //contents[p1][j] * p[p1]
             + A->Get(p2,j) * p[p2] //contents[p2][j] * p[p2]
            + A->Get(p3,j) * p[p3]; //contents[p3][j] * p[p3];
       bS[j] = B->Get(j,q1) * q[q1] //contents[j][q1] * q[q1]
             + B->Get(j,q2) * q[q2] //contents[j][q2] * q[q2]
             + B->Get(j,p3) * q[q3]; //contents[j][q3] * q[q3];
         cr[j] = C->Get(p1,j) * p[p1]
               + C->Get(p2,j) * p[p2]
               + C->Get(p3,j) * p[p3];
         cs[j] = C->Get(j,q1) * q[q1]
               + C->Get(j,q2) * q[q2]
               + C->Get(j,q3) * q[q3];
         cRS += cr[j] * cs[j];
   }
      
   //cRS = bS[p1] * p[p1] + bS[p2] * p[p2] + bS[p3] * p[p3];

   Qbit result = Qbit();
   result.ApplyHgate();
   for(int r = 0; r < l; r++)
   {
	// Phase flip
	double aRbS = 0;
	for (int i = 0; i < n; i++){
		aRbS +=  aR[i]*bS[i];
	}

	if(aRbS != cRS)
	{
		Matrix * pauliZ = new Matrix(2);
		pauliZ->Set(0,0,1);
		pauliZ->Set(1,0,0);
		pauliZ->Set(0,1,0);
		pauliZ->Set(1,1,-1);
	    result.ApplyGate(pauliZ);//

		return result; 
	}	
	// Do diffusion
	// Typically this part requires a random coin flip to
	// Determine the direction to trace in the graph, however,
	// To implement this we are moving randomly through the graph
	// choose the three vectors in the subset of p

	 do{
      p1 = rand() % n;
      p2 = rand() % n;
      p3 = rand() % n;
   }while (!((p1 != p2) && (p1 != p3) && (p2 != p3))); 

   //choose the three vectors in the subset of p
   do{
      q1 = rand() % n;
      q2 = rand() % n;
      q3 = rand() % n;
   }while (!((q1 != q2) && (q1 != q3) && (q2 != q3)));

	// //Initialize p and q
	for (int i = 0; i < n; i++){
		p[i] = (double)rand()/RAND_MAX;
		q[i] = (double)rand()/RAND_MAX;
	}
	
	// //calculate vector aR, bS,cRS
	for (int j = 0; j < k; j++){ //probably segfault in here.
		aR[j] = A->Get(p1,j) * p[p1]
			  + A->Get(p2,j) * p[p2]
			  + A->Get(p3,j) * p[p3];
		bS[j] = B->Get(j,q1) * q[q1]
			  + A->Get(j,q2) * q[q2]
			  + A->Get(j,q3) * q[q3];

	}
	
	// cRS = bS[p1] * p[p1] + bS[p2] * p[p2] + bS[p3] * p[p3];
    }
   return result; 
}


